library("affy")
library("geneplotter")
library("SpikeInSubset")
library("limma")
data(spikein95)

pd <- data.frame(population=c(1,1,1,2,2,2),
replicate=c(1,2,3,1,2,3))
rownames(pd) <- sampleNames(spikein95)
vl <- list(population="1 is control, 2 is treatment",
           replicate="arbitrary numbering") 
phenoData(spikein95) <- new("phenoData",pData=pd,varLabels=vl)


##some probe level stuff
pms <- pm(spikein95)
mms <- mm(spikein95)
pns <- probeNames(spikein95)
smoothScatter(log2(mms[,1]),log2(pms[,1]))
abline(0,1)

Index <- split(seq(along=pns),pns)
mypar(2,3)
par(ask=TRUE)
for(i in 1:10){
  for(j in 1:6){
    matplot(log2(t(rbind(pms[Index[[i]],j],mms[Index[[i]],j]))),type="b",pch=c("P","M"),main=names(Index)[i])
  }
}
###summarize
eset <- rma(spikein95)
design <- model.matrix(~factor(eset$population))
fit <- lmFit(eset, design)
ebayes <- eBayes(fit)
tab <- topTable(ebayes,coef=2,adjust="fdr",n=10)


genenames <- as.character(tab$ID)



library("hgu95av2")
hgu95av2()
gn <- mget(genenames, hgu95av2GENENAME)
sapply(gn, function(x) if(!is.na(x)) strwrap(x) else x)



ls(hgu95av2GENENAME)[1:10]


library("GO")
library("annotate")
gos <- mget(genenames, hgu95av2GO)
gos[[1]][1:2]
goids <- sapply(gos, names)
goids[[1]]
goterms <- sapply(goids, function(x) mget(as.character(x), GOTERM))
sapply(goterms[[1]], Term)



library("XML")
absts<-pm.getabst(genenames,"hgu95av2")



absts[[1]][[4]]



###################################################
##  commented out pm.titles since it behaves weird. 
## Revert back to it when it is fixed?
##titl <- pm.titles(absts[2])
titl <- sapply(absts[[2]], articleTitle)
strwrap(titl, simplify=FALSE)


pro.res <- sapply(absts, function(x) pm.abstGrep("[Pp]rotein", x))
pro.res[[2]]



pmAbst2HTML(absts[[2]],filename="pm.html")



ll <- getLL(genenames,"hgu95av2")
ug <- unlist(lookUp(genenames, "hgu95av2", "UNIGENE"))
sym <- getSYMBOL(genenames,"hgu95av2")



tab <- data.frame(sym,tab[,-1])
htmlpage(list(ll, ug),filename="report.html",title="HTML report",othernames=tab,
         table.head=c("Locus ID", "UniGene ID", colnames(tab)),table.center=TRUE, 
          repository=list("ll","ug"))



library("annaffy")
atab <- aafTableAnn( genenames, "hgu95av2", aaf.handler() )
saveHTML(atab, file="report2.html")



library("GOstats")

index <- sample(1:12000, 100)
gn <- ls(hgu95av2GENENAME)[index]
lls <- unique(getLL(gn, "hgu95av2"))
hyp <- GOHyperG(lls, "hgu95av2", what = "MF")
index2 <- hyp$pvalues < 0.05
wh <- mget(names(hyp$pvalues)[index2], GOTERM)
terms <- sapply(wh, Term)
strwrap(terms[1:10], simplify = FALSE)
hyp$pvalues[index2][1:10]


